% Scripts to demonstrate the effect of red noise on tidal line estimations 
% 1) create a red spectrum with of frequency slope -2 2) add tidal lines and 2) invert it using LS to get the
% tidal amplitudes & phase back. 

%% Background spectra of slope -2
%define coefficients

A = [1 -.999]; % slope  of -2
numsamples = 1024*100;
nfft = 1024;
fs = 1/60; % sampling every 5 minutes

% create background red signal. Adjust the amplitudes to control the slope
% effect
background = filter(1,A,randn(numsamples,1));

% Add daily variations. 
Amp = [1 0.75 0.5 0.1];
Frq = [1/24 1/12 1/8;]./3600; % 24 hur, 12 hr and 8 hour
t = (1:numsamples)/fs;
y =     Amp(1)*sin(2*pi*t*Frq(1)+pi/4) ...
    +   Amp(2)*sin(2*pi*t*Frq(2)+pi/4) ...
    +   Amp(3)*sin(2*pi*t*Frq(3)+pi/4) ...
    +   normrnd(0,Amp(4),1,length(t));

% amplitude modulation y = x.*cos(2*pi*fc*t)
modfreq = 2/(365.25 * 24 * 3600 );  % semi-annual modulation
y = y.* abs(Amp(1)*cos(2*pi*t*modfreq));


% Add them to create synthetic signals
sdata =  background + y';

% now plot the power spectra
[Pyy,F] = pwelch(sdata, hanning(nfft),1,nfft,fs);
loglog(1./(F*3600),abs(Pyy),'b', 'LineWidth', 2);


%% Recover the amplitudes
periods = [24.0000 12.0000 8.0000 ].*3600;

% periods = [7.0000 17.0000 23.0000].*3600; some random periods

%make a model


m=[cos(2*pi*t(:)/(periods(1)) ) sin(2*pi*t(:)/(periods(1)) )];

for j=2:length(periods)
   m=[m cos(2*pi*t(:)/(periods(j)) ) sin(2*pi*t(:)/(periods(j)) )]; 
end

% now invert with the whole synthetic signals

 [b,stats] = robustfit(m,sdata); % data with "red" noise
 [a, stats] = robustfit(m,y); % data without "red" noise
 
 a_amps = abs(complex(a(2:2:end),a(3:2:end)));
 b_amps = abs(complex(b(2:2:end),b(3:2:end)));
 
 

 fprintf('Expected White Red\n');
 for i = 1:3,
 fprintf('%6.4f %6.4f %6.4f\n', Amp(i),a_amps(i),b_amps(i));
 end;

 
 %% pre whitening 
 
% pre whitening is a process to make the signal "white". Typically the geomagnetic 
% data are "red" meaning, the power increases with period. We should make
% the data as white as possble before inverting the data to get tidal
% amplitudes. This is done by 1) fitting the tidal lines to the data 2)
% compute residuals and 3) find a filter that makes the residul white 4)
% apply that filter to the original data 5) repeating the steps so that the
% residuals become complete white.
 


 % Get aplitude and phase 
P = angle(complex(b(2:2:end),b(3:2:end)));
A = abs(complex(b(2:2:end),b(3:2:end)));

for i = 1:length(P),
   yy(i,:) =   A(i) * sin(2*pi*t./periods(i) + P(i));
end;


IR = cra([sdata-sum(yy)' sum(yy)'], 2500);
C = conv(sdata, IR);

[d, stats] = robustfit(m,C'); % robustfit
d_amps = abs(complex(d(2:2:end),d(3:2:end)));


 fprintf('Expected Red PW\n');
 for i = 1:3,
 fprintf('%6.4f %6.4f %6.4f\n', Amp(i),b_amps(i),d_amps(i));
 end;


